  function [hnew,ctrim,hchol]=hessian_three( h )  ; 
% function [hnew,ctrim,hchol]=hessian_three( h )  ; 
%
% Input
% ------
% H n x n matrix need not be pd 
% 
% Use [Hc,d]=chol(H) and extract sequentially matrices of maximum dimension
% that are positive definite 
% For entries corresponding to d, set to abs val of the diagonal entries in the original 
% matrix 
% 
% Output 
% ------
% hnew      PD matrix 
% ctrim     Count of trimmings needed to obtain PD matrix 
% hchol     hchol'*hchol = hnew 
% 
n=size(h);
if n(1)~=n(2) 
    error('Input h must be square matrix')
end 
n=n(1); 

hnew=[]; 
htrim = h; 

ctrim = 0; 

istop = 0; 
while istop < 1 
    
    [hc,d] = chol( htrim ); 
    
    if d == 0 
        hnew=blkdiag( hnew, hc ) 
        istop = 1 
        
    else 
        ctrim = ctrim + 1 ;
        nt = size( htrim , 1 );
        htemp = zeros( d ); 
        
        if d > 2 
            htemp(1:d-1,1:d-1) = hc ; 
        end 
        
        htemp( d, d) = sqrt( abs(htrim(d,d)) ); 
        
        if d < nt 
            htrim=htrim(d+1:end,d+1:end); 
        else 
            istop = 1; 
        end 
        
        hnew = blkdiag( hnew , htemp ); 
        clear htemp nt; 
        
    end 
    
end 
if size( hnew, 1) ~= n 
    error('Code produces matrix of wrong size') 
end 
hnew = hnew'*hnew; 
[hchol,p] = chol( hnew ); 
if p~=0 
    error('function hessian_tree.m still does not produce pd matrix') 
end